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Conversion of Osculating Orbital Elements to Mean Orbital Elements 
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Abstract 


Orbit determination and ephemeris generation or prediction over relatively long elapsed times can be 
accomplished with mean elements. The most simple and efficient method for orbit determination, which is also 
known as epoch point conversion, performs the conversion of osculating elements to mean elements by iterative 
procedures. Previous epoch point conversion methods are restricted to shorter elapsed times with linear 
convergence. The new method presented in this paper calculates an analytic initial guess of the unknown mean 
elements from a first order theory of secular perturbations and computes a transition matrix with accurate 
numerical partials. It thereby eliminates the problem of an inaccurate initial guess and an identity transition 
matrix employed by previous methods. With a good initial guess of the unknown mean elements and an 
accurate transition matrix, converting osculating elements to mean elements can be accomplished over long 
elapsed times with quadratic convergence. 


Basic Concepts 

This paper presents new methods to solve the following problems: 

• A user’s propagator requires a mean orbital element set (e.g., NORAD-North American Aerospace 
Defense Command two-card element set or any other mean element set) as input but the element set is not 
available. 

• For ground based radar acquisition and space sensor surveillance, a single osculating state vector of an 
object at the current time is available, but the mean elements corresponding to an epoch a few days, weeks 
or months earlier are not. If the set of mean elements at an epoch can be computed, then the object can be 
identified with respect to a known catalog (e g. NORAD element set for Resident Space Objects). The 
mean elements at an epoch are needed to efficiently provide radar or sensor pointing commands. 

The osculating orbital elements represent, in a general sense, the true position and velocity vectors of a satellite, 
but are poorly behaved over time as a basis for prediction. The mean orbital elements do not represent the true 
position and velocity vectors of a satellite, but are well behaved over time. An orbit described by a set of mean 
orbital elements is said to be perturbed or non-Keplerian. 

By way of notation all vectors are in bold unless specified. For the sake of simplicity, the osculating elements 
and the mean elements are respectively denoted as: 

y(t) = [a e i Q © M] t and y(t) = [a e i fl a m[ 

where a is the semimajor axis, e is the eccentricity, i is the inclination, Q is the longitude of the ascending 
node, © is the argument of perigee and M is the mean anomaly. The six elements of y(t) or y(t) can be chosen 
in a variety of ways and the classical orbital elements are chosen to enhance theoretical understanding. The 
singularities of small eccentricity, small inclination or critical inclination do not exist in the conversion of 
osculating elements to mean elements, but the reverse is not true. The transformations between classical orbital 
elements, y(t), and predicted or true position and velocity vectors, r(t) and v(t), are simple. The transformation 
of the mean elements, y(t), to the mean position and velocity vectors, r(t) and v(t), is straightforward, but 
the reverse is difficult. By way of definition, a transformation between vectors is an instantaneous conversion. 
A conversion between vectors involves an elapsed time greater than or equal to zero. 
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0 A 

The ballistic coefficient of a satellite is defined as B = — , where Cj® is the zero drag coefficient, A is the 

2m 

reference area of the satellite and m is the mass of the satellite. The problem of converting osculating orbital 
elements at a final time, t, to mean orbital elements at an initial time, tg, can be stated as 

Given: B, t 0 , t and y(t) 

Find: y( t Q ) 

Depending on the perigee altitude of the satellite, previous methods presented by References 1 and 2 that use 
iterative procedures, are restricted to a short elapsed time,(t- t 0 ). Reference 3, which uses a combination of 
Kozai’s and Izsak’s theories, calculates the difference between the osculating elements and mean elements with 
zero elapsed time. The difference, y( t 0 ) - y( t 0 ), is the very small variation of orbital elements due to short- 

periodic perturbations. The method presented in this paper allows the elapsed time to be extended over much 
longer intervals. Only the case of non-negative elapsed time, (t— t 0 ) £ 0, is formulated, however, the method 

also holds for negative elapsed time as well. 

To understand our approach to solve the above problem, a background on some basic concepts is required. The 
state prediction problem for a satellite orbiting about a central body such as the Earth is to find the position and 
velocity vectors, r(t) and v(t), at time t that satisfy the vector equation of motion 


subject to the given initial conditions r = r 0 and v = v 0 at time t = t 0 . In Equation (1), p is the gravitational 
constant and a d is the total disturbed acceleration vector due to disturbed gravity, atmospheric drag, lunisolar 
gravitational attractions, solar radiation pressure, tidal friction, n-body gravitational attractions and thrust. The 
first term on the right-hand-side of Equation (1) is the acceleration vector due to central gravity. 

If a d is zero, then Equation (1) can be solved analytically by one of Kepler’s methods and the osculating 

elements (a, e, i, Q, <b) describing the size, shape and orientation of the satellite orbit remain constant in time. 
The osculating mean anomaly, M, defines the angular position of the satellite in its orbit with respect to time. 
For a typical non-thrusting, near-Earth satellite between the theoretical atmospheric altitude of 91 km ^(300,000 

ft) and a 12-hour orbit altitude of approximately 20,000 km, the distuibed acceleration vector, a d> is due 

mainly to the Earth zonal gravitational harmonics, J 2 , J 3 and J 4 > and atmospheric drag. Lunisolar 
gravitational attractions and solar radiation pressure, whose effects may be formulated similar to atmospheric 
drag, are neglected in this study. The orbital elements affected by the distuibed acceleration vector, a*, are 
defined as the osculating elements, y(t). [Strictly speaking the osculating elements are affected by a d .] The 
disturbed accelerations of a* cause secular, short-periodic and long-periodic variations in the classical oibital 

elements (a, e, i, £2, <d, M). Short periods are on the order of time of one satellite passage around the Earth. 
Long periods are on the order of time of one complete perigee passage around the Earth. If the periodic effects 
are removed, then the new orbital elements are defined as mean elements, y(t) . That is, the mean elements are 

affected only by secular perturbations. 

Periodic variations occur in all osculating elements and are induced by all zonal gravitational harmonics. 
However, the variations in osculating elements induced by periodic perturbations are much smaller than those 
ind uced by secu l a r pertuibations as the elapsed time increases. The variations in osculating elements induced 
by secular pertuibations are constant or non-periodic, and only even zonal gravitational harmonics and 
atmospheric drag give rise to secular effects. Atmospheric drag can be a significant part of the secular 
pertuibations if the perigee altitude of a typical satellite orbit is less than 500 km. In summary: 

• Secular pertuibations » Periodic pertuibations (for long elapsed time) 

• Secular perturbations = Distuibed acceleration due to gravity ( J 2 and J 4 ) 

and Drag (perigee altitude < 500 km) 
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Figures 1 to 6 show the variations of the osculating elements due to the disturbed acceleration vector, a d , for a 

non-thrusting Low Earth Orbit satellite with a perigee altitude of approximately 200 km. Osculating elements 
are computed by four methods: numerical integration, two NORAD propagators (SGP and SGP4) and a first 
order theory of secular perturbations. The osculating elements predicted by the first order theory of secular 
perturbations are depicted by the thick solid line. As shown in Figures 1 to 3, the first order secular effects on 
the osculating elements a, e and i are almost negligible and are induced only by drag, and therefore the 

averaged time derivatives — , — and — are almost zero (the thick solid lines are almost horizontal). 

dt av dt av dt av 

The mean elements, y ( t 0 ) , are the initial values of osculating elements at time t 0 = 0 on the thick solid line. 
The solutions computed by numerical integration, SGP and SGP4, are depicted respectively by the thin solid 
line, the triangle and the square. 


As shown in Figures 4 to 6, the secular effects on the osculating elements Q, co and M are significant and are 


induced by both gravity and drag. The averaged time derivatives 


an 

dt 


dco 

"dt" 


and 


dM 

dt 


are almost constant 


to first order even if the time dependent contributions due to drag are included. At 200 km altitude, the 
disturbed acceleration due to J 2 is at least one or two orders of magnitude greater than that due to J 4 It is well 
known in general perturbations theory that the osculating elements (a, e, i, Q, co) are “slow” variables and that 

dMl 


M is a “fast” variable. This implies that 


dt 


dn 

is much greater than — ■ 
dt | 

dM 


and — | 
dt 


If the elapsed time. 


(t- t 0 ), is long, then the averaged time derivative due to drag, — — , must be included for a typical Low 

dt 


drag 


Earth Orbit satellite with a perigee altitude of 200 km even though References 4 to 6 and many excellent 
textbooks have recommended otherwise. 


Since the satellite orbit of this example is almost circular, the first order estimates of the osculating eccentricity 
and argument of perigee are not equal to their averaged values during the two periods of the satellite orbit as 
shown in Figures 2 and 5. However, the first order estimates of the other osculating elements are very close to 
the averaged values of the osculating elements as predicted by the NORAD propagator SGP. 

If a NORAD propagator is used, then the mean mean motion, n , which replaces a of y(t), must be carefully 
computed. The osculating mean motion, n, is determined from the equation: p = n a . If the elapsed time is 
short (less than a day), the errors arising from interchanging the osculating mean motion and mean mean 
motion are negligible. If the elapsed time is long (on the order of days), the mean mean motion, if, must be 
used in evaluating the averaged time derivatives (except for initialization). It should be clear that p * n a . 


From general perturbations theoiy, the osculating elements and the mean elements have a first order secular 
relationship given by 


dy i 

y(t) = y(t 0 ) + -g 


0 -‘ 0 ) 


where the averaged time derivatives valid between t 0 and t, are defined as 
dy 


dt 


1 (ju 1 2 b dy 

= — | — dM = — f (1-ecosE) — dE 

271 o dt 2 7t o dt 


( 2 ) 


( 3 ) 


dy , 

and E is the eccentric anomaly. The instantaneous time derivative, — , can be obtained from Lagrange s 

dt 

dy 

planetary equations for secular perturbations. Analytic solutions of — 


, which are not difficult to obtain, are 


relatively accurate except for Low Earth Orbit satellites. If the elapsed time is long and the disturbed 
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acceleration due to drag is significant (e g. for a Low Earth Orbit satellite), then the analytic averaged time 
derivatives of the semimajor axis and eccentricity can be an order of magnitude in error with respect to their 
numerically integrated values. This can happen even if the elapsed time is less than one period for some 
satellite orbits. 

From Reference 7, the averaged time derivatives due to the disturbed gravity (J 2 and J 4 ) for secular 
perturbations are given by 


da 

cU 

de 

dt 

di 

dt 


= 0 


gravity 


= o 


gravity 


= 0 


gravity 


dQ 


3ncos i r 2 

b-fi 

^12-80sin 2 i^-e 2 ^4 + 15sin 2 ij| 

dt 

gravity 

2p 2 

5J 4 re 

(4-7sin 2 iV2+3e 2 )] 




16p 2 

[\ A /J j 

do 



35, » r -U-Wil 


dt 

gravity 

4p 2 l 

J 



( 4 ) 

( 5 ) 

( 6 ) 

( 7 ) 


+ — jlOsin 2 i^76-89sin 2 ij+e 2 ^56-36sin 2 i- 45 sin 4 i^J 


> nJ 2 r 4 
384p 
15n J 4 r 4 


( 8 ) 


32p 4 


dM 

dt 


= n + 


^16-62 sin 2 i + 49sin 4 ij + e 2 ^18-63sin 2 i+-^^sin 4 ij 
3n 0 J 2 r 2 (l-e 2 ) 


gravity 


4p 2 




1/2 


-|2-3sin 2 ij 

sin 2 i c (lOO — 13 1 sin 2 ij+c 2 ^20- 98sin 2 i+67sin 4 ijj 
e 2 ^8 -40 sin 2 i +35 sin 4 ijj 


( 9 ) 


96p 4 

45nJ 4 r 4 (l-e 2 ) 

128p 4 

where p = a(l-e 2 ) . In satellite state prediction, an accurate solution must always be computed numerically. 

Therefore, the averaged time derivatives due to drag should be integrated numerically for the semimajor axis 
and eccentricity. From References 8, the averaged time derivatives due to drag for secular pertuibations are 
given by 

dal 


dt 


= -2Bna 2 


drag 


1 71 o 

— J pQi(l + e cosEX 
71 o 




ecosE) 

ccosE) 


dE 


= - 2Bn p 
d[ 

dt drag 


-J pQiQ 2 M^ 

*0 V0- e 


cosE) 


cosE) 


dE 


= 0 


( 10 ) 

( 11 ) 

( 12 ) 
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(13) 


<m 


dt 

drag 

d© 


dt 

drag 

dM 


dt 

drag 


= 0 
= 0 


= i„_ v , .21* 

2 ' w 4a dt 


drag 


(t-t 0 ) 


where the rotating Earth factors are: 


Q 0 = -^2® 


Qi = 1 - Qo 


i Vl-e 2 


(1 — e cosE) 
(1 + e cosE) 


Qn — — — 7 

0 2 = cosE - (1 - e cosE) (2 cosE - e - e cos E) 

2(1 -e 2 ) 


(14) 

(15) 


and co e is the constant scalar Earth rotational rate and p is the intantaneous density at a reference altitude with 
respect to the mean anomaly E. 

The drag-induced averaged time derivatives of inclination, longitude of the ascending node and argument of 
perigee are normally much smaller than those of their gravity counterparts, and therefore are neglected. The 
averaged second time derivatives of the semimajor axis and eccentricity are at least three orders of magnitudes 
less than the averaged time derivatives and therefore are also neglected. 

Since the disturbed accelerations are additive, the averaged time derivatives of the osculating elements y(t) due 
to secular perturbations are the sum of the gravity and drag components. Substituting the averaged time 
derivatives of Equations (4) through (15) into Equation (2) and rearranging, a good initial guess of the 
unknown mean dements y ( t Q ) at time t 0 is given by 


' a( V ' 


a(t) 

e(t 0 ) 


e(t) 

i(‘o) 


i(0 

«<t 0 ) 


fi(t) 

®( t 0 ) 


©(0 

M( t 0 ) 


M(t) 


da 

dT 


( 


dQ| 

dt 


dt 

0 

>(t- W 


(t-t 0 ) 

drag 

de 

- <*- *0> 
drag 


dco 

( "dT 


gravity 


)(t-t 0 ) 


gravity 


( 


dMj 

dt 


gravity 


dMj 

dt 


drag 


)(t-t 0 ) 


(16) 


The mean dements a , e , i and n in the averaged time derivatives of Equations (4) through (15) are initially 
replaced by the respective osculating elements to start the initial guess process of Equation (16). Fortunately, 
the variations of these mean elements are normally slow. Now, given B, t 0 , t and y(t), the right-hand-side of 
Equation (16) can be computed to give a good initial guess of the unknown mean dements y( t 0 ) at time t Q . 
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Methods for Converting Mean Elements to Osculating Elements 

A method of general perturbations for satellite theory seeks the solution of Equation (1) by series expansion and 
term-by-term analytic integration of the disturbed acceleration. General perturbations methods circumvent 
numerical integration but must initiate with mean elements. The NORAD element set propagators, which are 
general perturbations methods, start with a given mean vector, y( t 0 ) = (n 0 ,e 0 ,i 0 ,f^,co 0 ,M 0 ) at a given 

epoch time, t Q . The given mean mean motion, n 0 , can be converted to the mean semimajor axis, a 0 , as 
described in Reference 9. Reference 10 documented five NORAD models for propagating Resident Space 
Objects (RSO) and satellites around the Earth. The Simplified General Perturbations (SGP) model, which 
contains most of the first order gravitational terms as described in References 11 and 12, computes the drag 
terms as linear functions of time. The SGP4 model, which uses the gravitational model of References 13 and 
14, calculates the drag terms by a power density function of the atmosphere. The SGP4 model is customarily 
used for near-Earth satellites. A space object is classified by NORAD as near-Earth if its period is less than 225 
minutes; it is classified as deep-space otherwise. The SDP4 model, which includes the gravitational terms due 
to third-body effects of the Sun and Moon and the Earth sectorial and tesseral harmonics, is an extension of 
SGP4 for deep-space objects. The SGP8 model is an extension of SGP4 with the same gravitational and drag 
models, but predicts state vectors more accurately especially when the satellite altitude is under 200 km. The 
SDP8 model is an extension of SGP8 and SDP4 for deep-space objects. The two higher order propagators, 
SGP8 and SDP8 were briefly considered as replacements for the SGP4 and SDP4, but the increased 
computational time and only slight improvement in state prediction accuracy have discouraged the changeover 
until the present time. 

New and improved versions of SGP4 and SDP4 may be obtained directly from NORAD. However, using a 
UNIX workstation which is linked to Internet, a version of the five NORAD propagators (SGP, SGP4, SDP4, 
SGP8 and SDP8) can be downloaded from a computer at the Air Force Institute of Technology. The 
comprehensive instructions of Reference 10 and all the necessary algorithms downloaded from Internet do not 
guarantee that the reader can use the NORAD propagators immediately. Reference 9 discusses the problems 
and solutions. 


The NORAD propagators were developed from the satellite theories of Kozai and Brouwer to propagate mean 
elements to osculating elements. Other mean element to osculating element propagators such as those 
described in References 16 and 17 will not be considered since their improvements in computational speed, 
memory storage and state vector accuracy are insignificant for our purpose. A method of converting osculating 
elements to mean elements requires a “forward” propagator which converts mean elements to osculating 
elements. Figures 1 to 6 show that the predicted osculating elements propagated by SGP are closer to those 
predicted by the first order theory of secular perturbations (the thick solid line) for most satellite orbits. If the 


forward propagator is SGP, then — = J 

2 4a dt 1 


drag 


, it 5n da 

and — = 

6 12a dt 


are required and can be 


'drag 


computed from Equation (10). 


If an accurate — 
6 


is required, then 

dr 


needs to be computed. If the 


drag 


forward propagator is SGP4 or SDP4, then B*gp 4 = 6378137.0 K p, B, where 2£ K £ 1 is a constant 


related to the density model used to calculate p 0 at the altitude of 120 km. The units of p 0 and B are 
respectively kg/ m 3 and m 2 / kg . As a concrete example, a SGP or NORAD propagator has been chosen as 
the forward propagator in this paper even though the choice is aibitraiy. 


The NORAD element set propagators all start with a given mean vector, y( t 0 ), and their output is the 
predicted position and velocity vectors, r(t) and v(t), which can then be transformed to the osculating elements, 
y(t). Osculating elements are the ones that are usually available and the reconstruction of mean elements must 
begin with osculating elements. 
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Previous Methods for Converting Osculating Elements to Mean Elements (Ref. 1 to 3) 


A method for converting osculating elements to mean elements which uses a transition matrix 


g y(0 

# y( t 0 ) 


( 17 ) 


for the equation 

8y(t) = T8y(t 0 ) (18) 

is known as a method of differential corrections or a transition matrix algorithm in applied optimal control 
theory. A transition matrix algorithm usually begins with a guessed nominal y( t Q ) at time t 0 , and then 

propagates forward to a nominal y(t) at time t, which in turn gives 8y( t). Therefore the differential corrections 
at time t 0 is 


Sy( to) = T- 1 8y( t) (19) 

using Equation (18). In practice, accurate transition matrices are computed numerically. Chapter 7 of 
Reference 18 provides the algorithm to numerically compute a transition matrix such as that of Equation (17). 
Reference 19 describes the technique of accurate numerical partials. 


Reference 1 suggests an iterative procedure by using the given osculating elements y(t) as the initial guess of 
the mean elements y( t 0 ) and assuming an identity transition matrix. That is, the averaged time derivatives of 
Equation (16) are zero for any elapsed time and Equation (19) is reduced to 

5y(t 0 ) = 8y(t) (20) 

since T = I. If the elapsed time is less than one day, this iterative procedure may convert y(t) to y( t 0 ) for some 
Low Earth Orbits. This method converges linearly at best. If the elapsed time is greater than a day, this 
method fails for most Low Earth Orbits. The cause of failure is a combination of: 

• The neglected drag terms for semimajor axis, eccentricity and mean anomaly are not small. 

• The identity matrix is a poor approximation to the transition matrix and Equation (20) is not valid. 


Recalling that the problem is to find y( t 0 ) given B, t 0 , t and y(t). The iterative procedures of Reference 1 
may be summarized as follows: 

1. Let the initial guess of the mean elements at the given time t 0 be the same as the given osculating 
elements, y(t). That is 

ykOo) = y(‘) 

where k is the iteration number (k = 0 at this point). 

2. Propagate forward (using a SGP propagator) from y k ( t 0 ) to x k ( t) and then transform x k ( t) to 
y k ( t) . The difference in osculating elements at time t is 

SyO) = y(t) - ykO) = SyCV 

using Equation (20). 

3. Compute the new guess of the mean elements at time t 0 as 

yk+iOo) = yk( l 0 ) + g y( V 

For 8y = | y( t) - y k ( t) | > 10“ 10 , then procedure 2 is repeated with the new guess 
y k+ i( t Q ); otherwise the desired mean elements at time t Q is 

y( l 0 ) = yk+iOo) 


References 2 and 3 are transformations between osculating and mean elements with t = t Q . Reference 2 is an 

iterative method that uses the Frazer elements (position and velocity vectors and their variations), includes the 
short- and long-periodic perturbations. Reference 3, which uses a combination of Kozai’s and Izsak’s theories, 
calculates the difference between the osculating elements and mean elements only due to short-periodic 
perturbations. The methods of these two references are only instantaneous conversions or transformations, and 
therefore are not described in this paper. 
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A New Method for Converting Osculating Elements to Mean Elements 

This method uses a good initial guess derived from a first order theory of secular perturbations and a transition 
matrix computed by accurate numerical partials. As shown in Figures 1 to 6, the osculating elements computed 
from a first order theory of secular perturbations (the thick solid straight line) behave linearly with respect to 
time and are close to the average values of the SGP solutions. This implies that the initial guess of the mean 
elements computed from Equation (16) will be close to the desired mean elements, y(t 0 ), at time t 0 . A 
traditional transition matrix algorithm requires 7 forward propagations (1 nominal and 6 neighboring 
trajectories). The stepsize hj is usually set to 10" 6 of the i ^ mean element; this may be too small for 
eccentricity but too large for the semimajor axis. The transition matrix algorithm based on accurate numerical 
partials requires 25 forward propagations (1 nominal and 24 neighboring trajectories). The iterative procedures 
of this method may be summarized as follows: 


1 . 


2 . 


3. 

4. 

5. 


Compute the averaged time derivatives of Equation (4) to (15). Let the given osculating elements, y(t), 
be Y, then the initial estimate of the mean semimajor axis, mean eccentricity, mean inclination and 


mean mean 

motion are given by: 



a 

= a(t) - 

dt 

(t- 

drag 

■to) 

e 

= e(t) - 

dt 

(t- 

drag 

-to) 

i 

= i(0 




n 

= n(t) 





where the empirical constant X is 0.5 for Low Earth Orbits and zero otherwise. 

Initialize — and — for the forward propagator SGP and the initial guess of the mean mean anomaly, 
2 6 


M, 


by 


numerically 


da, 

integrating — | 


of Equation 


( 10 ), 


then 


■drag 


n 

2 


, h 5n da 

— and— = 

Hdrag 6 12a dt 


drag 


. If the forward propagator is SGP4 or SDP4, then 


3n da 
4 a dt 
is required. 

Compute the nominal mean elements y*( t 0 ) guess from Equation (16) and then propagate forward to 
t giving the nominal osculating elements y*( t) 

Compute the nominal differential correction of osculating elements at t as 
8y(t) = Y - y*(t) 

r d y(t) 


Compute the transition matrix, T = 


by accurate numerical partials. 


d y(t 0 ) 

Propagate forward by a SGP propagator 6 times in the neighborhood of the nominal trajectory (step 
3) using 

y(t 0 ) = y*(t 0 ) + 8y(t 0 ) 


For the 1 th neighboring trajectory ( i = 1, 2, .. , 6 ), four “neighboring” neighboring trajectories are 
computed from 4 stepsizes of y, -y, -y- and --y- . The corresponding 4 osculating elements 

computed by a SGP propagator at time t are yj, yj and y^ . The partial derivatives of the i* 
column of T is given by 
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6. 


£ y(t) = (y3-y4)-p 3 (yi~y2) 
s yi(‘o) ph,(i- p 2 ) 


The 6x6 transition matrix can be approximated as 


T = 

0 y(t) 


d y(t) d y(t) 

0 y(t) 


. 0 y(t 0 ) 


0 yi(t 0 ) 0 j^Oo) 

0 ynflo) . 


In computing T, 24 propagation by the forward propagator SGP are required. A good choice of p is 
1/2 and that of the stepsize hj is io - 4 of each mean element. 


Update the nominal y* ( t 0 ) at t 0 as 

y*O 0 )| = y*U 0 )| 


lold 


+ T" 1 5y( t) 


Step 3 to 6 are repeated until the magnitude of 5y( t) is reduced to an acceptably small value (10 12 ). 


Examples 


Two examples are given to illustrate the performance of the new method for satellites at a Low Earth Orbit and 
a High Earth Orbit. Using a first order theory of secular perturbations, the initial guess of the slow variables, 
(a,e,i , Q , co ) , can be estimated in the vicinity of their unknown values at time t 0 . If the initial guess of the 

fast variable, M , can be predicted to within approximately 30 degrees of its unknown value, then convergence 
is fast. This is not a problem for almost any satellite orbit if the elapsed time is less than one day. 


The variations of the osculating mean anomaly, M, at time t for the example Low Earth Orbit satellite are 
shown in Figures 6. The initial guess of M at time t 0 is related to M at time t by Equation (16). For short 

elapsed times, the 30 degree requirement can be satisfied easily. The primary reason is that the unknown mean 
semimajor axis, a , has changed only slightly in a short elapsed time, and as a consequence the effects on the 


terms due to 


dM 

dt 


gravity 


and 


dM 

dt 


drag 


are small. The right-hand-side of Equation (16) can then be computed 


quite accurately giving a good initial guess of M . 


In the following examples, first we assumed to know the mean elements, y( t 0 ), at time 1$ , and then we used a 

forward propagator (SGP4 for Example 1 and SDP4 for Example 2) to get the osculating elements, y(t), at time 
t. In what follows, we discard the mean elements, y( t 0 ), and no knowledge of the mean elements at time 

will be used. The problem is: 

Given B, t 0 , t and y(t); retrieve y( t 0 ). 


Example 1 

The mean elements, y( t 0 ) , at time t 0 are taken out of the SGP4 example of Reference 10 with minor 
adjustments for double precision computation. Also the ballistic coefficient, B, is replaced by that of the 
LANDSAT-D. The osculating elements, y(t), at time t are computed by using SGP4 with 
from Equation (18). 

C A 

Given: C do = 2.0, A = 12.2778 meter 2 , m = 1710.0 kg, B = = 0.00718 

2m 

t 0 = 0.0, t = 0, 1, 5 days 


B*gp 4 computed 

meter 2 

kg 
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'6641.774062 

.0096661858 

72.85385095 

115.9622955 

59.40458042 

103.8371428 


y(l day) = 


y(5 days) = 


6641.774062 6633.640850 6589.666059 

.0096661858 .0083375605 .0053578232 

72.85385095 72.85513567 j 72.84433117 

y(0 day) = n5 9622 955 ’ y(lday)_ 113.4116703 ’ y(5days)_ 103.0554340 
59.40458042 56.77943987 54.07490762 

103.8371428 130.9131346 349.1054119 

Following the procedures 1 to 3 of the new method, the initial guesses of the nominal mean elements, y*( tg) , 
at time t 0 for the cases of t = 0, 1, 5 days are computed as: 

f 664 L77401 f6638.463ll f664 1.92571 


'6589.666059' 

.0053578232 

72.84433117 

103.0554340 

54.07490762 

349.1054119 



y (V = 

1 0 day 

Comparison of results: 


Days from t Q 


664L7740 

.00966618 

72.853850 

115.96229 

59.404580 

103.83714 


y <M. 


.00869400 

72.855135 

115.97053 

59.232901 

105.09619 


y*( V 


6641.9257 

.00834486 

72.844331 

115.93401 

66.308201 

138.72967 


Method of Reference 1 (Walter) 

New Method (Der&Danchick) 

# of iterations 
required: Nj 

# of SGP4 calls 7 

N, 

# of iterations 
required: N 2 

# of SGP4 calls 25 

n 2 

12 

84 

3 

75 

40 

280 

3 

75 



One reason to change the ballistic coefficient from that of Reference 10 is to investigate how close the new 
method works during the satellite orbital decay. Using the lifetime equation of Reference 4, this satellite has a 
lifetime of approximately 9000 minutes or a little over six days from time t 0 . Numerical integration shows 

that the lifetime is approximately 10 days from time t 0 . After 5 days, the perigee altitude of the satellite is 

close to 150 km and the state vectors predicted by the SGP4 propagator become inaccurate. Nevertheless, the 
new method converges very close to y( t 0 ) in the last few days before satellite re-entry. 


Example 2 

The mean elements, y( t 0 ), at time t 0 are constructed from a Molniya orbit with perigee altitude of 500 km 

and apogee altitude of 40000 km. Critical inclination (63.4 degrees) is chosen for this HEO to demostrate that 
there is no singularity at any inclination for the conversion from osculating elements to mean elements. The 
ballistic coefficient, B, which is not important for this HEO, is chosen to be the same as that of example 1. The 
osculating elements, y(t), at time t are computed by using the SDP4 propagator. The effect of drag is 
negligible on this HEO satellite. 
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Given: C do = 2.0, A = 12.2778 meter 2 , m = 1710.0 kg, B = 

t 0 = 0.0, t= 10, 100, 200 days 

'26626.701651 [26630.49652' 

.7415398328 .7364330943 

x 63.45549845 63.36886054 

y( 10 days) = , y(100 days) = 

118.5145116 * y ’ 105.0926980 

.0695130482 .4435084366 

137.6237041 77.71324416 


= 0.00718 


kg 


y(100 days) = 


26630.49652 

.7364330943 

63.36886054 

105.0926980 

.4435084366 

77.71324416 


y(200 days) = 


26626.65222 

.7323070455 

63.72013696 

90.19640476 

.7794822322 

26.78924405 


Following the procedures 1 to 3 of the new method, the initial guesses of the nominal mean elements, y*( t^ , 
at time tg for the cases of t= 10, 100, 200 days are computed as: 


y*(t 0 ) 


26626.798' 
.74154077 
63.455498 
todays 119.93263 
.07822485 
144.12848 


y*( V 


26630.570 
.73643493 
63.368860 
too days 119.43259 
.43299544 
158.12109 


y (*o) 


26626.686 
.73230790 
63.720136 
200 days 117.80345 
1.5247843 
156.05599 


Comparison of results: 


Days from t Q 
(days) 


10 


100 


200 



Method of Reference 1 (Walter) 

# of iterations 
required: Nj 

# of SDP4 calls 7 

N, 

13 

91 

20 

140 


New Method (Der&Danchick) 


# of iterations 
required: N 2 


3 


# of SDP4 calls 25 
N, 



Conclusions 

• If the constants given in the data block of the five NORAD propagators are defined in double precision, 
then the osculating state vectors can be predicted much more accurately especially for satellite orbits 
without the influence of atmospheric drag. With this simple modification, the five NORAD propagators 
can be used as forward propagators for the conversion of osculating elements to mean elements. 

• The method presented by this paper achieves quadratic convergence due to accurate numerical partials. It 
has uniformly good performance over the three test cases, succeeds where the Reference 1 method failed, 
and is generally more computationally efficient. 
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• For Low Earth Orbits (the example described by Figures 1 to 6), the initial guess of the fast changing mean 
mean anomaly, M , at time t 0 must include the term computed from the averaged time derivative due to 

drag. In this case, the conversion of osculating elements to mean elements can be extended from an 
elapsed time of one day to the last few days before satellite re-entry 

• For High Earth Orbits, the initial guess of the fast changing mean mean anomaly, M , at time t 0 can be 

predicted accurately even for long elapsed times. In this case, the elapsed time for the conversion of 
osculating elements to mean elements can be extended to months. 

• For Geosynchronous Earth Orbits (results not included in this paper), the initial guess of the fast changing 
mean mean anomaly, M , at time t 0 can also be predicted accurately even for long elapsed times. In this 

case, the elapsed time for the conversion of osculating elements to mean elements can be extended to years. 

• The determination of mean orbital elements can be radically streamlined with the new method and made 
applicable for mean elements orbit determination with observations made by heterogeneous sensor types 
over long spans of elapsed time. 
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Figure 1. Effects of gravity and drag on semimajor axis for 2 Low Earth orbits 
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Figure 2. Effects of gravity and drag on eccentricity for 2 Low Earth Orbits 





Figure 3. 
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Figure 4. Effects of gravity and drag on 


33 < 




Mean anomaly (degrees) Argument of perigee (degrees) 



Figure 5. Effects of gravity and drag on argument of perigee for 2 Low Earth Orbits 



Figure 6, Effects of gravity and drag on mean anomaly for 2 Low Earth Orbits 
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